package scu.maqiang.fes;

import scu.maqiang.mesh.Mesh;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SIMatrix;
import scu.maqiang.numeric.SRMatrix;

import java.util.Arrays;

public class FES2T3Mix_RT0P0 extends FES {

    public FES2T3Mix_RT0P0(Mesh2T3 mesh) {
        super(mesh, 1);
        edges = mesh.renumberingEdges();
        //System.out.println(edges);
        nEdge = edges.getNNZ();
        int[] edgeI = new int[nEdge];
        int[] edgeJ = new int[nEdge];
        int[] edgeNum = new int[nEdge];
        edges.getTriValues(edgeI, edgeJ, edgeNum);
        int nt = mesh.getNt();
        ndof = nEdge + nt;
        edgeElements = new int[nt][];
        for(int i = 0; i < nt; i++) {
            int[] ele = mesh.getElement(i);
            int[] edgeEle = new int[3];
            edgeEle[0] = Math.max(edges.getElement(ele[0], ele[1]), edges.getElement(ele[1], ele[0])) - 1;
            edgeEle[1] = Math.max(edges.getElement(ele[1], ele[2]), edges.getElement(ele[2], ele[1])) - 1;
            edgeEle[2] = Math.max(edges.getElement(ele[2], ele[0]), edges.getElement(ele[0], ele[2])) - 1;
            edgeElements[i] = edgeEle;
        }
        //System.out.println(Arrays.deepToString(edgeElements));
    }
    int[][] edgeElements;
    int nEdge;
    SIMatrix edges;

    public int getNumberOfEdges() {
        return nEdge;
    }
    public int getEdgeNumber(int i, int j) {
        return Math.max(edges.getElement(i, j), edges.getElement(j, i)) - 1;
    }

    public void extractSolution(double[] sol, double[] ux, double[] uy, double[] u) {
        double[] allArea = new double[mesh.getNv()];
        for(int i = 0; i < mesh.getNt(); i++) {
            int[] ele = mesh.getElement(i);
            int[] sign = new int[3];
            sign[0] = ele[0] < ele[1]? 1:-1;
            sign[1] = ele[1] < ele[2]? 1:-1;
            sign[2] = ele[2] < ele[0]? 1:-1;
            mesh.getCoordPerElement(ele, coordRegion);
            double[] lenEdge = new double[3];
            double area = CG.area2DT3(coordRegion);
            double[] x = MVO.col(coordRegion, 0);
            double[] y = MVO.col(coordRegion, 1);
            CG.edge2DT3(coordRegion, lenEdge);
            int[] dof = edgeElements[i];
            double[] localSol = new double[3];
            for(int j = 0; j < 3; j++) {
                localSol[j] = sol[dof[j]];
            }
            double p1x = (sign[0] * localSol[0] * lenEdge[0] * (x[0] - x[2]) + sign[2] * localSol[2] * lenEdge[2] * (x[0] - x[1])) / (2 * area);
            double p1y = (sign[0] * localSol[0] * lenEdge[0] * (y[0] - y[2]) + sign[2] * localSol[2] * lenEdge[2] * (y[0] - y[1])) / (2 * area);

            double p2x = (sign[0] * localSol[0] * lenEdge[0] * (x[1] - x[2]) + sign[1] * localSol[1] * lenEdge[1] * (x[1] - x[0])) / (2 * area);
            double p2y = (sign[0] * localSol[0] * lenEdge[0] * (y[1] - y[2]) + sign[1] * localSol[1] * lenEdge[1] * (y[1] - y[0])) / (2 * area);

            double p3x = (sign[1] * localSol[1] * lenEdge[1] * (x[2] - x[0]) + sign[2] * localSol[2] * lenEdge[2] * (x[2] - x[1])) / (2 * area);
            double p3y = (sign[1] * localSol[1] * lenEdge[1] * (y[2] - y[0]) + sign[2] * localSol[2] * lenEdge[2] * (y[2] - y[1])) / (2 * area);

            double[] px = {p1x, p2x, p3x};
            double[] py = {p1y, p2y, p3y};
            for(int j = 0; j < 3; j++) {
                ux[ele[j]] += px[j] * area;
                uy[ele[j]] += py[j] * area;
                u[ele[j]] += sol[nEdge + i] * area;
                allArea[ele[j]] += area;
            }
        }
        for(int i = 0; i < mesh.getNv(); i++) {
            ux[i] /= allArea[i];
            uy[i] /= allArea[i];
            u[i] /= allArea[i];
        }
    }

    public void assembleMixedStiff(double constCoef, BVPType tp, SRMatrix A) {
        int nt = mesh.getNt();
        double[][] B = new double[3][3];
        double[] C = new double[3];
        for(int i = 0; i < nt; i++) {
            int[] ele = mesh.getElement(i);
//            System.out.println("========================");
//            System.out.println(Arrays.toString(ele));
            int[] sign = new int[3];
            sign[0] = ele[0] < ele[1]? 1:-1;
            sign[1] = ele[1] < ele[2]? 1:-1;
            sign[2] = ele[2] < ele[0]? 1:-1;
//            System.out.println(Arrays.toString(sign));
//            System.out.println(MVO.toString(edgeElements));
            int[] dof = edgeElements[i];
//            System.out.println(Arrays.toString(dof));
            //计算得到刚度矩阵中的B与C

            mesh.getCoordPerElement(ele, coordRegion);
//            System.out.println(MVO.toString(coordRegion));
            double[] lenEdge = new double[3];
            double area = CG.area2DT3(coordRegion);
            CG.edge2DT3(coordRegion, lenEdge);
//            System.out.println(Arrays.toString(lenEdge));
            double[] x = MVO.col(coordRegion, 0);
            double[] y = MVO.col(coordRegion, 1);
            double mx = MVO.average(x);
            double my = MVO.average(y);
            double x2y2 = x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[0] * x[1] + x[1] * x[2] + x[2] * x[0];
            x2y2 += y[0] * y[0] + y[1] * y[1] + y[2] * y[2] + y[0] * y[1] + y[1] * y[2] + y[2] * y[0];
            x2y2 /= 6.0;
//            System.out.println("x2y2 = " + x2y2);
            for(int ii = 0; ii < 3; ii++) {
                double[] pi = coordRegion[(ii + 2) % 3];
                double edgeI = lenEdge[ii];
                for(int jj = 0; jj < 3; jj++) {
                    double[] pj = coordRegion[(jj + 2) % 3];
                    double temp = x2y2 - (mx * (pi[0] + pj[0]) + my * (pi[1] + pj[1])) + pi[0] * pj[0] + pi[1] * pj[1];;
                    B[ii][jj] = 0.25 * edgeI * lenEdge[jj] / area * temp / constCoef;
                }
                C[ii] = sign[ii] * edgeI;
            }
//            System.out.println(MVO.toString(B));
//            System.out.println(Arrays.toString(C));
//            System.exit(1);
            for(int ii = 0; ii < 3; ii++) {
                //C[ii] *= sign[ii];
                for(int jj = 0; jj < 3; jj++) {
                    B[ii][jj] *= sign[ii] * sign[jj];
                }
            }

//            System.out.println(MVO.toString(B));
//            if (i == 10) {
//                System.exit(1);
//            }
            A.assemble(dof, B);
            A.assemble(i + nEdge, dof, C);
            A.assemble(dof, i + nEdge, C);
        }
    }

    public void assembleMixSource(double const_dirichletBC, double const_Source, double[] RHS, int... bdLabel) {
        int nt = mesh.getNt();
        double[] RHS_dirichletB = new double[2];
        double RHS_source;
        for(int i = 0; i < mesh.getNb(); i++) {
            int label = mesh.getBoundaryLabel(i);
            for(int j = 0; j < bdLabel.length; j++) {
                if (label == bdLabel[j]) {
                    int[] ele = mesh.getBoundary(i);
                    mesh.getCoordPerBoundary(ele, coordBoundary);
                    int dof = getEdgeNumber(ele[0], ele[1]);
                    double len = CG.length2DL2(coordBoundary);
                    RHS[dof] = RHS[dof] + len * const_dirichletBC;
                    break;
                }
            }
        }

        for (int i = 0; i < mesh.getNt(); i++) {
            int[] ele = mesh.getElement(i);
            mesh.getCoordPerElement(ele, coordRegion);
            double area = CG.area2DT3(coordRegion);
            int dof = i + nEdge;
            RHS[dof] = RHS[dof] - const_Source * area;
        }

    }
}
